# coding: utf-8
import itertools
import numpy as np
from mammoth.core import Bond, Angle, Dihedral, Improper, ClassicalForceField
__all__ = ["QM2FF"]
[docs]class QM2FF:
"""
A class that implements the Seminario method to get a force field
(1) Bautista, E. J.; Seminario, J. M. Int. J. Quantum Chem. 2008, 108 (1), 180–188.
(2) Seminario, J. M. Int. J. Quantum Chem. 1996, 60 (7), 1271–1277.
TODO: update the doc with new references when the methods will be improved
**General attributes**
.. attribute :: molecule
A molecule object which contains at least the hessian matrix.
.. attribute :: natom
Number of atom in the molecule.
.. attribute :: bond_list
List of tuples of two atom indexes which define a bond according to
the bond orders or the connectivity.
*Warning:* depending on the bond search strategy, the bond_list and bonds
are not consistents because the bond_list only contains bond defined by
the bond_orders or the connectivity.
.. attribute :: bonds
A list of Bond objects which describe a bond with atom indexes, bond length,
force constants and atomic species.
.. attribute :: angle_list
List of tuples of 3 atom indexes which define an angle.
.. attribute :: angles
A list of Angle objects which describe an angle with atom indexes,
angle value in degree, force constants and atomic species.
.. attribute :: dihedral_list
List of tuples of 4 atom indexes which define a dihedral angle.
.. attribute :: dihedrals
A list of Dihedral objects which describe a dihedral angle with atom
indexes, dihedral angle value in degree, force constants and atomic species.
The periodicity can be computed if neighbors are provided.
.. attribute :: improper_list
List of tuples of 4 atom indexes which define an improper torsion.
.. attribute :: impropers
A list of Improper objects which describe an improper torsionwith atom
indexes, torsion angle value in degree, force constants and atomic species.
.. attribute :: connectivity
The connectivity of the molecule as a list of atom index pairs. The connectivity
depends on the bond search strategy.
**Parameters for the calculation**
.. attribute :: cutoffb
Cutoff distance in angstrom for in order to build the bond list.
.. attribute :: bond_search
Bond search strategy :
* hessian (default): bond exists if the force constant is real
* bond_order: bond exists if it is in the molecule.bond_orders dictionary
* hybrid: bond exists if the force constant is real or if it is in the
molecule.bond_orders dictionary
* connectivity: consider only bonds defined by the connectivity. In that
case, the connectivity must be provided.
.. attribute :: imag_th
Threshold in order to determine if an imaginary part is zero.
.. attribute :: dihedral_type
Dihedral type is harmonic or parabolic depending on the functional form
of the potential energy. Default is harmonic.
*Note about complex eigenvalues :* If the sub hessian matrix of used to
compute the force constant of each internal coordinates cannot be diagonnalized
in R, numpy return the complex eigenvectors and complex eigenvalues. Thus
the test concerning the imaginary part of the eigenvalues is relevant even if
the hessian matrix is not explicitely a complex matrix.
"""
def __init__(self, molecule, connectivity=None, cutoffb=2.5,
bond_search="hessian", imag_th=1e-3, dihedral_type="harmonic",
improper_type="all", cutoff_sp2=10.):
"""
Set up QM2FF object to compute a force field
Args:
molecule (Molecule): A molecule object with hessian matrix and bond orders
cutoffb (float): cutoff distance for bond search
bond_serch (str): Define how bonds are identified. Authorized values
are "hessian" (default), "bond_order" or "hybrid"
imag_th (float): threshold under which imaginary part is zero (default 0.001)
connectivity (list): if bond_search is 'connectivity' you have to provide
the molecular connectivity as in a PDB file:
[(0, 1), (1, 2, 3), (2, 5), ...]
where the first index is bonded to the following
indexes.
dihedral_type (str): 'harmonic' or 'parabolic'
improper_type (str): 'all' or 'sp2'
cutoff_sp2 (float): if improper_type is sp2, cutoff ...
"""
self.molecule = molecule
self.natom = len(self.molecule)
# parameters
self.cutoffb = float(cutoffb)
self.bond_search = bond_search.lower()
self.imag_th = float(imag_th)
self.dihedral_type = dihedral_type
self.connectivity = connectivity
self.improper_type = improper_type
self.cutoff_sp2 = float(cutoff_sp2)
# TODO: These threshold are internal attribute. Make it accessible from outside ?
self._align_threshold = 1.e-2 # threshold for dot products
# check parameters
self._check_params()
# set up all attributes
self.bond_list = None
self.bonds = None
self.angle_list = None
self.angles = None
self.dihedral_list = None
self.dihedrals = None
self.improper_list = None
self.impropers = None
def _check_params(self):
""" check parameters for the FF calculations """
# Define here all tests to do when the QM2FF object is instanciated
# this will allow to check if all parameters are right and present
# before the calculation is run.
if self.bond_search not in ("hessian", "hybrid", "bond_order", "connectivity"):
raise ValueError("The bond search strategy is wrong: %s" % self.bond_search)
if self.bond_search == "connectivity" and not self.connectivity:
raise ValueError("The bond search strategy is 'connectivity' " +
"and the connectivity was not provided.")
if self.dihedral_type not in ("parabolic", "harmonic"):
raise ValueError("The dihedral type is wrong: %s" % self.dihedral_type)
try:
if not self.molecule.bond_orders:
raise ValueError("A bond orders dict is needed.")
except AttributeError:
print("Use the core.forcefield.Molecule class with a filled bond_orders dict.")
raise AttributeError("Molecule object has no attribute bond_orders")
def _get_connectivity(self):
"""
Return the connectivity as a list of atom index pairs. If the connectivity
was not defined, the returned connectivity is the one obtained from
the chosen bond search strategy. In consequence, the connectivity may
differ from the bond orders.
"""
if self._connectivity:
return self._connectivity
else:
if not self.bonds:
self.get_bonds()
return set([b.atoms for b in self.bonds])
def _set_connectivity(self, connectivity):
""" set connectivity and do some check """
if not connectivity:
self._connectivity = None
else:
# force a list of tuples corresponding to atom index pairs
bonds = []
for pair in connectivity:
if len(pair) != 2:
print(pair)
raise ValueError("The connectivity is not a list of pairs")
elif not all(isinstance(iat, int) for iat in pair):
print(pair)
raise ValueError("The connectivity is not a list of integer pairs")
else:
bonds.append(tuple(pair))
bonds = set(bonds)
# check for double entry
remove = list()
for iat, jat in bonds:
if (jat, iat) in bonds:
remove.append((jat, iat))
[bonds.remove(bond) for bond in remove]
self._connectivity = bonds
connectivity = property(_get_connectivity, _set_connectivity)
[docs] def set_PDB_connectivity(self, connectivity):
"""
Convert the connectivity initially in a format such as PDB in a list of
atom index pairs and store the resulting list in a set object.
WARNING: atom indexes in PDB usually start at 1. Thus atom indexes are
translated before being store in the connectivity.
Example of connectivity in PDB format:
CONNECT 1 3
CONNECT 2 3
CONNECT 3 1 2 4
CONNECT 4 3 5 6 7
CONNECT 5 4
CONNECT 6 4
CONNECT 7 4
A list such as the following is expected:
[[1, 3], [2, 3], [3, 1, 2, 4], [4, 3, 5, 6, 7], [5, 4], [6, 4], [7, 4]]
"""
bonds = list()
for sublist in connectivity:
iat = sublist[0]
bonds += [(iat - 1, jat - 1) for jat in sublist[1:]]
bonds = set(bonds)
# check for double entry
remove = list()
for iat, jat in bonds:
if (jat, iat) in bonds and (iat, jat) not in remove:
remove.append((jat, iat))
[bonds.remove(bond) for bond in remove]
self._connectivity = bonds
[docs] def get_PDB_connectivity(self):
"""
Return the connectivity in a PDB format. The connectivity depends on the
chosen bond search strategy. The output string is such as :
CONNECT 1 3
CONNECT 2 3
CONNECT 3 2 4 1
CONNECT 4 3 7 5 6
CONNECT 5 4
CONNECT 6 4
CONNECT 7 4
"""
lines = ""
for iat in range(self.natom):
begin = "CONNECT%5d" % (iat + 1)
line = ""
for i, j in self.connectivity:
if i == iat:
line += "%5d" % (j + 1)
elif j == iat:
line += "%5d" % (i + 1)
if line != "":
lines += begin + line + "\n"
return lines
[docs] def get_forcefield(self):
""" Compute all parameters and return a ClassicalFrorceField object """
# compute parameters :
self.get_bonds()
self.get_angles()
if self.dihedral_type == "harmonic":
self.get_dihedrals_harmonic()
elif self.dihedral_type == "parabolic":
self.get_dihedrals_parabolic()
self.get_impropers()
cFF = ClassicalForceField(self.molecule, self.bonds, self.angles,
self.dihedrals, self.impropers,
self.molecule.at_charges)
return cFF
[docs] def get_bonds(self):
"""
Build the list of bonds according to the bond search strategy and compute
the associated force constants. Bond search strategy are :
* hessian: bond exists if the force constant is real
* bond_order: bond exists if it is in the molecule.bond_orders dictionary
* hybrid: bond exists if the force constant is real or if it is in the
molecule.bond_orders dictionary
* connectivity: only consider bond defined by the connectivity
The distance cutoff between atom is defined by self.cutoffb.
"""
self.bonds = []
self.bond_list = []
for iat, jat in itertools.combinations(range(self.natom), 2):
# check if the bond order is not null
in_bond_orders = False
if ((iat, jat) in self.molecule.bond_orders or
(jat, iat) in self.molecule.bond_orders):
in_bond_orders = True
# check if the bond is defined in the connectivity
in_connectivity = False
if self._connectivity:
if (iat, jat) in self._connectivity or (jat, iat) in self._connectivity:
in_connectivity = True
if self.bond_search == "bond_order" and not in_bond_orders:
# in bond_order case, consider only bond in the bond_order dict
continue
if self.bond_search == "connectivity" and not in_connectivity:
# in connectivity case, consider only bond in the connectivity list
continue
# compute distance and bond vector
v_ij = self.molecule[jat].coords - self.molecule[iat].coords
r_ij = np.linalg.norm(v_ij)
v_ijN = v_ij / r_ij
if r_ij > self.cutoffb:
if self.bond_search == "hessian":
# in hessian case do not consider bond longer than cutoffb
continue
elif self.bond_search == "hybrid" and not in_bond_orders:
# in hybrid case do not consider bond longer than cutoffb
# if the bond order is null
continue
# compute force constant
H = self.molecule.hessian[3 * jat:3 * jat + 3, 3 * iat:3 * iat + 3]
eigenValues, eigenVectors = np.linalg.eig(H)
eigvN = [eigenVectors[i] / (np.linalg.norm(eigenVectors[i])) for i in range(3)]
Kb = sum([abs(np.dot(v_ijN, eigvN[i])) * eigenValues[i] for i in range(3)])
if not self._check_exist(eigenValues) or Kb < 0.:
if self.bond_search == "hessian":
# invalid bond in hessian case
continue
elif self.bond_search == "hybrid" and not in_bond_orders:
# invalid bond in hybrid case if the bond order is null
continue
elif self.bond_search in ("bond_order", "connectivity"):
print("WARNING: wrong force constant for bond ", (iat, jat))
print("Kb = ", Kb, type(Kb))
abond = Bond(atoms=(iat, jat), value=r_ij, frc_cste=Kb.real,
elements=(self.molecule[iat], self.molecule[jat]))
self.bonds.append(abond)
# in the bond list, only bonds defined by the bond order or the connectivity
# are stored. This list will be used to identify angles, dihedrals and
# impropers.
if self.bond_search == "connectivity":
self.bond_list = set(self.connectivity[:])
else:
self.bond_list = set(self.molecule.bond_orders.keys())
[docs] def get_angle_list(self):
"""
Build the list of possible angles according to the bond order list or
the connectivity. In consequence the angle list does not depend on the
bond search strategy.
Angles are store in the `angle_list` attribute.
"""
# build the bond list if needed
if not self.bond_list:
self.get_bonds()
self.angle_list = list()
for iat, jat, kat in itertools.permutations(range(self.natom), 3):
if iat < kat:
if (iat, jat) in self.bond_list or (jat, iat) in self.bond_list:
if (jat, kat) in self.bond_list or (kat, jat) in self.bond_list:
self.angle_list.append((iat, jat, kat))
self.angle_list = set(self.angle_list)
return self.angle_list
[docs] def get_dihedral_list(self):
"""
Build the list of possible dihedral angles according to the bond order
list or the connectivity. In consequence the dihedral angle list does
not depend on the bond search strategy.
Dihedral are store in the `dihedral_list` attribute.
"""
# build the bond list if needed
if not self.bond_list:
self.get_bonds()
paras_list = list()
for iat, jat, kat, lat in itertools.permutations(range(self.natom), 4):
if iat < lat:
if (iat, jat) in self.bond_list or (jat, iat) in self.bond_list:
if (jat, kat) in self.bond_list or (kat, jat) in self.bond_list:
if (kat, lat) in self.bond_list or (lat, kat) in self.bond_list:
paras_list.append((iat, jat, kat, lat))
self.dihedral_list = set(paras_list)
return self.dihedral_list
[docs] def get_improper_list(self):
"""
Build the list of all possible improper torsions according to the bond
order list or the connectivity. In consequence the improper list does not
depend on the bond search strategy.
Improper torsions are store in the `improper_list` attribute.
"""
# build the bond list if needed
if not self.bond_list:
self.get_bonds()
# first, build the whole list
paras_list = list()
for iat, jat, kat, lat in itertools.permutations(range(self.natom), 4):
if jat < kat and kat < lat:
if (iat, jat) in self.bond_list or (jat, iat) in self.bond_list:
if (iat, kat) in self.bond_list or (kat, iat) in self.bond_list:
if (iat, lat) in self.bond_list or (lat, iat) in self.bond_list:
paras_list.append((iat, jat, kat, lat))
# second, do some check and filtering
remove = list()
add = list()
for improper in paras_list:
# chek if j, k, l are aligned
if self._check_aligned(*improper[1:]):
print("Improper ", improper)
print("WARNING: j, k and l are aligned. I do not consider that improper")
remove.append(improper)
continue
# check if i, j, k are aligned
if self._check_aligned(*improper[0:3]):
# issue if i, j, k are aligned
# if r_mid < 1e-1 A, i is between j and k => reverse k and l
remove.append(improper)
iat, jat, kat, lat = improper
improper = (iat, jat, lat, kat)
add.append(improper)
if self.improper_type == "sp2":
# in that case, if i is far away the plane defined by atoms
# (j, k, l), the improper is not considered.
v_ij = self.molecule[improper[0]].coords - self.molecule[improper[1]].coords
v_jk = self.molecule[improper[1]].coords - self.molecule[improper[2]].coords
v_kl = self.molecule[improper[2]].coords - self.molecule[improper[3]].coords
uN = np.cross(v_kl, v_jk)
uN /= np.linalg.norm(uN)
d = np.abs(np.dot(uN, v_ij))
if d > self.cutoff_sp2:
# i, j, k, l are not in the same plane => remove that improper
remove.append(improper)
# remove/add impropers according to the above filtering
[paras_list.remove(improper) for improper in remove]
[paras_list.append(improper) for improper in add]
# store the list
self.improper_list = set(paras_list)
return self.improper_list
[docs] def get_angles(self):
"""
Compute the force constant associated to angles store in the angle_list
attribute and sotre the results in a list of Angle objects.
"""
if not self.angle_list:
self.get_angle_list()
self.angles = list()
for angle in self.angle_list:
vv = [self.molecule[angle[1]].coords - self.molecule[angle[0]].coords,
self.molecule[angle[1]].coords - self.molecule[angle[2]].coords]
rr = [np.linalg.norm(v) for v in vv]
vvn = [v / r for v, r in zip(vv, rr)]
H = [self.molecule.hessian[3 * angle[1]: 3 * angle[1] + 3,
3 * angle[i]: 3 * angle[i] + 3] for i in [0, 2]]
eigenValues, eigenVectors = np.linalg.eig(H)
if not self._check_exist(eigenValues):
print("angle ", angle)
print("WARNING: the force constant of that angle presents " +
"an imaginary part. You should check the geometry.")
angle_val = np.degrees(np.arccos((np.dot(vv[0], vv[1])) / (rr[0] * rr[1])))
eigvN = [[eigenVectors[i][j] /
np.linalg.norm(eigenVectors[i][j]) for j in range(3)] for i in range(2)]
un = np.cross(vvn[1], vvn[0]) / (np.linalg.norm(np.cross(vvn[1], vvn[0])))
up = [np.cross(un, vvn[0]), np.cross(vvn[1], un)]
Kap = [(rr[i] ** 2) * sum([abs(np.dot(up[i], eigvN[i][j])) * eigenValues[i][j]
for j in range(3)]) for i in range(2)]
Ka = 1 / ((1 / Kap[0]) + (1 / Kap[1]))
if Ka > 0:
an_angle = Angle(atoms=angle, value=angle_val, frc_cste=Ka.real,
elements=tuple(self.molecule[i] for i in angle))
self.angles.append(an_angle)
[docs] def get_impropers(self):
"""
Compute the force constant associated to improper torsions store in the
improper_list attribute and sotre the results in a list of Improper objects.
"""
if not self.improper_list:
self.get_improper_list()
self.impropers = list()
remove = list()
for improper in self.improper_list:
# 0 1 2 3 4 5
#[v_ij, v_ik, v_il, v_jk, v_jl, v_kl]
vv = [self.molecule[improper[i]].coords - self.molecule[improper[j]].coords
for i in range(3) for j in range(4) if i < j]
rr = [np.linalg.norm(v) for v in vv]
vvn = [v / r for v, r in zip(vv, rr)]
H = [self.molecule.hessian[3 * improper[0]:3 * improper[0] + 3,
3 * improper[i + 1]:3 * improper[i + 1] + 3]
for i in range(3)]
eigenValues, eigenVectors = np.linalg.eig(H)
if not self._check_exist(eigenValues):
print("Improper ", improper)
print("WARNING: the force constant of that improper presents " +
"an imaginary part. You should check the geometry.")
eigvN = [[eigenVectors[i][j] / (np.linalg.norm(eigenVectors[i][j]))
for j in range(3)] for i in range(len(H))]
uN = [np.cross(vvn[3], vvn[0]), np.cross(vvn[5], vvn[3])]
uN = [u / np.linalg.norm(u) for u in uN]
Kan = sum([sum([eigenValues[x][i] * abs(np.dot(uN[1], eigvN[x][i]))
for i in range(3)]) for x in range(3)])
omega_r = np.arccos(np.dot(uN[1], uN[0])) # en radian
# np.acos(x) Return the arc cosine of x, in radians.
omega_d = np.degrees(omega_r) # en degré
# np.cos(x) Return the cosine of x radians
h_ijkl = np.dot(uN[1], vv[0]) # r_mid * np.cos(omega_r)
Ki = Kan * (h_ijkl**2)
# TODO: check here that h_ijkl is the right quantity
if Ki > 0:
an_improper = Improper(atoms=improper, value=omega_d, frc_cste=Ki.real,
elements=tuple(self.molecule[i] for i in improper))
self.impropers.append(an_improper)
else:
remove.append(improper)
# remove impropers to keep list consistents
[self.improper_list.remove(improper) for improper in remove]
[docs] def get_dihedrals_parabolic(self):
"""
Compute the force constant associated to dihedral angles considering a
parabolic expression of the potential energy along the coordinate as proposed
in the Seminario method (IJQC 30, 1996, 1271-1277).
V_phi = 1 / 2 k_phi (phi - phi_0)^2
All dihedral angles stored in the dihedral_list attribute are considered and
the force constants are store in a list of Dihedral objects.
"""
if not self.dihedral_list:
self.get_dihedral_list()
self.dihedrals = list()
remove = list()
for dihedral in self.dihedral_list:
# chek if i, j, k are aligned
if self._check_aligned(*dihedral[0:3]):
print("Dihedral ", dihedral)
print("WARNING: i, j, and k are aligned. I do not consider that dihedral")
remove.append(dihedral)
continue
# chek if j, k, l are aligned
if self._check_aligned(*dihedral[1:]):
print("Dihedral ", dihedral)
print("WARNING: j, k, and l are aligned. I do not consider that dihedral")
remove.append(dihedral)
continue
vv = [self.molecule[dihedral[i]].coords - self.molecule[dihedral[i + 1]].coords
for i in range(3)]
rr = [np.linalg.norm(v) for v in vv]
vvn = [v / r for v, r in zip(vv, rr)]
H = [self.molecule.hessian[3 * dihedral[i]:3 * dihedral[i] + 3,
3 * dihedral[i + 1]:3 * dihedral[i + 1] + 3]
for i in range(3)]
eigenValues, eigenVectors = np.linalg.eig(H)
if not self._check_exist(eigenValues):
print("dihedral ", dihedral)
print("WARNING: the force constant of that dihedral presents " +
"an imaginary part. You should check the geometry.")
eigvN = [[eigenVectors[i][j] / np.linalg.norm(eigenVectors[i][j])
for j in range(3)] for i in range(len(H))]
uN = [np.cross(vv[1], vv[0]) / np.linalg.norm(np.cross(vv[1], vv[0])),
np.cross(vv[2], vv[1]) / np.linalg.norm(np.cross(vv[2], vv[1]))]
dihedral_r = np.arccos(np.dot(uN[1], uN[0]))
cos = np.dot(uN[0], uN[1])
sin = np.dot(vvn[1], np.cross(uN[1], uN[0]))
tan = sin / cos
# TODO: je ne comprend pas la convention de signe
if sin < 0 and tan < 0:
dihedral_d = -1 * np.degrees(dihedral_r)
else:
dihedral_d = np.degrees(dihedral_r)
Kdp = [(rr[2 * i] ** 2) * (np.linalg.norm(np.cross(vvn[i], vvn[i + 1])) ** 2) *
sum([abs(np.dot(uN[i], eigvN[2 * i][j])) * eigenValues[2 * i][j]
for j in range(3)]) for i in range(2)]
Kd = 1 / (1 / Kdp[0] + 1 / Kdp[1])
if Kd > 0:
a_dihedral = Dihedral(atoms=dihedral, value=dihedral_d, frc_cste=Kd.real,
elements=tuple(self.molecule[i] for i in dihedral))
self.dihedrals.append(a_dihedral)
else:
print("dihedral ", dihedral)
print("Kd ", Kd)
print("WARNING: the force constant of that dihedral is negative")
remove.append(dihedral)
# remove dihedrals to keep list consistents
[self.dihedral_list.remove(dihedral) for dihedral in remove]
[docs] def get_dihedrals_harmonic(self):
"""
Compute the force constant associated to dihedral angles considering an
harmonic expression of the potential energy along the coordinate. The
force constant of the harmonic potential is computed from the force
constant of the parabolic potential. Assuming an harmonic potential
V_phi = K [1 + cos(n phi - d)]
The force constant K reads
K = K_p / n^2
where K_p is the parabolic force constant.
All dihedral angles stored in the dihedral_list attribute are considered and
the force constants are store in a list of Dihedral objects.
"""
# first, compute the parabolic dihedral force constant
if not self.dihedrals:
self.get_dihedrals_parabolic()
for dihedral in self.dihedrals:
jat, kat = dihedral.atoms[1:3]
nnj = sum(1 for bond in self.bond_list if jat in bond)
nnk = sum(1 for bond in self.bond_list if kat in bond)
dihedral.potential = "harmonic"
dihedral.compute_periodicity(nnj, nnk)
dihedral.frc_cste = dihedral.frc_cste / dihedral.periodicity**2
# TODO vérifier la périodicité
# TODO vérifier si la constante est bonne
def _check_exist(self, array):
"""
Return True if the imaginary part of all values in the array is lower
than the imag_th threshold
"""
return np.all(array.ravel().imag < self.imag_th)
def _check_aligned(self, i, j, k):
"""
Return true if the vectors ij and ik are aligned.
"""
vij = self.molecule[j].coords - self.molecule[i].coords
vij /= np.linalg.norm(vij)
vik = self.molecule[k].coords - self.molecule[i].coords
vik /= np.linalg.norm(vik)
aligned = False
if 1 - np.abs(np.dot(vij, vik)) < self._align_threshold:
aligned = True
return aligned